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1.    Introduction 

The  study  of  petroleum  reservoirs  is  characterized  by  strongly  nonlinear  equations, 
complex  physical  and  chemical  processes,  strong  spatial  variation  or  discontinuities  in  key 
reservoir  parameters,  uncertain  or  statistical  geological  data  and  unstable  fluid  regimes. 
Numerical  simulation  is  one  of  the  accepted  methods  for  the  study  of  petroleum  reservoirs 
and  improvements  in  numerical  methods  is  one  route  which  may  allow  progress  in  such  stu- 
dies. No  single  method  or  set  of  numerical  ideas  is  sufficient  at  the  present  time.  In  fact 
computational  simulation  is  used  for  many  distinct  length  scales,  and  to  suppress  or  represent 
accurately  a  wide  range  of  details  in  the  reservoir  and  fluid  description.  The  appropriate 
numerical  method  then  depends  on  the  level  of  description  required  and  the  purpose  of  the 
computation.  Similarly,  we  believe  that  a  variety  of  improved  methods  might  each  be  useful, 
possibly  in  distinct  contexts  or  facets  of  the  reservoir  simulation  problem.    In  the  same  vein. 
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we  believe  that  for  any  proposed  idea  or  method,  reservoir  parameters  and  computational 
problems  can  be  found  for  which  it  is  ill-suited. 

The  authors  and  coworkers  have  proposed  [2,10,11,13,12,8]  the  front  tracking  method 
as  useful  in  applications  to  petroleum  reservoir  simulation.  A  variety  of  tests  of  a  numerical 
analysis  nature  were  performed  for  the  method,  verifying  convergence  under  mesh  refine- 
ment and  absence  of  mesh  orientation  effects  [13].  The  ability  to  handle  complex  interface 
bifurcation  [8],  fingering  instabilities  [9,13]  and  polymer  injection  [3,4]  (as  an  example  of 
tertiary  oil  recovery)  indicates  a  level  of  robustness  in  this  method.  The  main  purpose  of  this 
paper  is  to  report  on  two  features  which  will  allow  further  series  of  tests  by  enabling  a  more 
realistic  description  of  reservoir  heterogeneities. 

The  reservoir  equations  and  the  numerical  method.  The  equations  governing  the  flow 
of  two  incompressible  phases  (oil  and  water)  in  a  porous  medium  can  be  approximated  by 

((t)  s),  +  V  •  V  f(.s)  -  0  ,  (1.1) 

V  =  -  K\  •  VP  ,  (1.2) 

V     V  =  V  •  K\  ■  VP  =  0  ,  (1.3) 

where  s  =  s(x,t)  is  the  fractional  volume  of  the  water  phase,  (J)(x)  is  the  porosity  of  the 
medium,  v(x.t)  is  the  total  fluid  (oil  plus  water)  velocity,  P(x,t)  is  the  pressure,  K(x)  is  the 
porous  medium  (rock)  permeability  tensor,  A.  =  \(s)  is  the  total  fluid  relative  transmissibil- 
ity  function,  and  f(s)  is  the  so-called  fractional  flow  function  relating  the  water  phase  velo- 
city to  the  total  fluid  velocity.  Equations  (1.1)  and  (1.3)  represent  conservation  of  the  fluids, 
(1.2)  is  Darcy's  law.  We  shall  consider  flow  in  two  spatial  dimensions.  In  order  to  concen- 
trate on  specific  physical  questions,  we  have  omitted  the  effects  of  gravity,  capillary  pressure 
(surface  tension),  variable  medium  depth,  and  flow  sources  from  (1.1)  through  (1.3). 
See  [24,25]    for  further  details. 

There  are  two  flow  regimes  associated  with  the  form  of  f(s).  The  case  of  a  fractional 
flow  function  linear  in  s  describes  miscible  flow;  a  non-linear  function  describes  immiscible 
flow.  For  miscible  flow  the  fluid  discontinuities  (shocks)  are  actually  contact  discontinuities 
and  the  shock  propagates  at  the  fluid  particle  velocity.  This  is  not  the  case  for  immiscible 
flow  and  in  that  case  the  fluid  particles  pass  through  the  shock  front.  As  a  consequence, 
though  the  hyperbolic  equation  (1.1)  has  a  simpler  wave  structure  for  miscible  flow,  it  is 
inherently  a  much  more  unstable  flow  regime  then  immiscible  flow.  In  the  linearized  (small 
perturbation)  regime,  the  stability  of  a  jump  discontinuity  for  (1.1)  -  (13)  is  shown  to  be 
determined  by  the  frontal  mobility  ratio 

Msb) 
m  =  -— —  (1.4) 

X(Sa) 

where  s^  (s^,)  is  the  state  on  the  ahead  (behind)  side  of  the  traveling  shock.  The  case  m  =  1 
is  the  limit  of  linear  stability,  with  m  >  1  corresponding  to  the  unstable  regime.  For  misci- 
ble flow  m  has  the  potential  of  becoming  infinite,  for  immiscible  flow  it  is  bounded  by  a 
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constant  value. 

The  analyses  in  this  paper  are  sparked  by  the  application  of  the  front  tracking  method  to 
reservoir  flow  problems.  Ihe  front  tracking  method  is  a  hybrid  which  combines  special 
adaptive  methods  for  the  enhanced  resolution  of  discontinuities  with  conventional  differenc- 
ing schemes  for  the  solution  in  the  region  between  discontinuities.  The  method  as  currently 
implemented  is  sequential;  the  pressure  (elliptic)  equation  (1.3)  and  the  conservation  (hyper- 
bolic) equation  (1.1)  are  solved  separately  every  time  step,  each  solution  involving  data  from 
the  previous  solution  of  the  other.  We  note  that  the  difference  in  characteristic  time  scales 
between  (1.1)  and  (1.3)  provides  sufficient  justification  for  this  splitting  of  the  system  (1.1), 
(1.3).  Employing  a  sequential  method  has  the  distinct  advantage  of  allowing  different  solu- 
tion techniques  to  be  used  for  each  equation.  The  pressure  equation  is  solved  by  the  method 
of  finite  elements,  which  is  well  suited  to  the  mathematical  character  of  elliptic  equations. 
Similarly  a  combination  of  finite  differences  and  analytical  Riemann  problem  solutions  is 
used  for  (1.1),  which  is  appropriate  in  view  of  its  hyperbolic  nature. 

The  presence  of  phase  or  other  types  of  discontinuities  in  the  physical  problem  gives 
rise  to  discontinuous  coefficients  in  the  elliptic  pressure  equation.  A  special  adaptive  grid, 
which  is  modified  at  each  time  step,  is  used  to  resolve  these  features  and  compute  the  flow 
field  accurately  [23]  .  Each  such  adaptive  grid,  which  we  shall  refer  to  as  the  'elliptic  grid', 
has  the  index  structure  of  a  regular  rectangular  Ng  x  Mg  grid  and  hence  the  numerical  solu- 
tion can  be  accelerated  by  fast  solution  techniques. 

In  the  hyperbolic  equation,  the  discontinuities  have  the  mathematical  structure  of  shock 
waves,  and  are  propagated  by  jump  relations  which  relate  the  shock  speed  to  the  magnitude 
of  certain  discontinuities  across  the  shock.  This  part  of  the  computation  can  be  viewed  as  a 
hybrid  of  finite  differences  and  of  the  method  of  characteristics  or  of  moving  point 
methods  [5].  For  a  discussion  of  front  tracking  in  greater  depth,  see  [2,13,14]  .  The  solu- 
tion of  (1.1)  in  the  regions  between  fronts  uses  a  finite  difference  scheme  with  respect  to  a 
regular  rectangular  Nh  x  Mj,  grid  which  covers  the  entire  computational  region,  and  which 
shall  be  referred  to  as  the  hyperbolic  grid. 

In  the  next  section,  we  describe  the  use  of  front  tracking  to  represent  geological  discon- 
tinuities such  as  layers  or  faults.  This  material  is  a  preliminary  report  on  a  portion  of  the 
Ph.D.  thesis  of  one  of  us  (M..VI.).  Next  we  report  on  the  effect  of  porosity  variation  on 
reservoir  fingering.  The  main  conclusion  is  that  porosity  is  less  significant  than  permeability 
as  a  cause  of  interface  instabilities.  A  final  section  contains  some  comments  on  .M.  Shearer's 
no-go  theorem. 
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2.   Geological  Layers 

2.1.    Introduction 

In  this  section  we  study  the  effects  of  a  discontinuous  permeability  tensor,  K,  brought 
about  by  the  presence  of  distinct  geological  layers.  We  concentrate  on  the  problem  of  a  sin- 
gle discontinuity  (a  zero  width  transition  zone)  separating  two  homogeneous  rock  layers  of 
different,  constant,  permeabilities.  We  assume  the  porosity  of  both  layers  is  the  same,  and 
constant,  and  scale  it  out  by  redefining  the  time  variable,  t.  We  are  interested  in  obtaining 
an  analytical  understanding  of  the  propagation  of  an  oil-water  phase  bank  (the  front)  through 
this  geological  discontinuity,  and  in  the  subsequent  development  of  a  numerical  algorithm  to 
embody  the  analytical  results  and  allow  calculations  through  media  of  greater  variation  in 
layering.  The  point  of  interaction  between  a  moving  front  and  the  stationary  rock  discon- 
tinuity will  be  referred  to  as  the  "node".  In  the  version  reported  on  here,  the  rock  discon- 
tinuity shall  be  assumed  to  be  either  horizontal  or  vertical. 

Such  an  interaction  falls  into  a  much  broader  class  of  interacting  discontinuities  of 
hyperbolic  systems.  In  such  interactions  one  is  interested  in  the  general  problems  of  bifurca- 
tion, deflection  and  evolution  of  the  intersection  point.  Glimm  and  Sharp  [6]  studied  this 
phase  bank  -  layer  discontinuity  interaction  as  an  example  of  elementary  waves  and  classified 
the  possible  exact  solutions  for  the  deflection  of  a  front  by  such  a  rock  discontinuity.  Assum- 
ing finite,  leading  order  data,  they  found  two  solutions.  The  first  consisted  of  a  one  parame- 
ter family  of  solutions  in  which  the  flow  is  parallel  to  the  front  and  the  node  remains  station- 
ary in  space.  The  second  is  a  solution  in  which  flow  is  normal  to  the  front,  and  hence  the 
node  propagates  in  space.  In  this  case,  the  angle  of  incidence  for  the  front  on  the  rock 
discontinuity  is  restricted  to  a  fixed  angle  given  in  terms  of  the  ratio  of  permeabilities  for  the 
two  geological  layers.  Both  of  these  solutions  are  too  restrictive  to  give  an  indication  of  how 
an  interaction  between  the  front  and  the  rock  discontinuity  develops,  though  they  give  possi- 
ble insight  into  steady  state  solutions.  A  third  solution  was  also  obtained  which  allowed  the 
shock  to  cross  the  layer  tangentially. 

This  problem  has  been  analyzed  with  greater  generality  [21].  The  theoretical  solution 
of  this  problem  for  general  angles  and  boundary  conditions  is  complicated  by  the  fact  that  the 
elliptic  equation  has  a  singularity  at  the  node  (the  velocity  is  (usually)  either  zero  or  infinity). 
Nevertheless  an  approximate  deflection  law  can  be  obtained  by  introducing  an  averaging 
length  scale  on  the  scale  of  the  ignored  physical  phenomena,  e.g.  capillarity.  The  resulting 
equations  will  be  discussed  on  another  occasion  [21).  In  the  next  section  we  describe  a  simple 
algorithm  that  will  give  an  approximate  method  for  deflection,  evolution,  and  lateral  bifurca- 
tion of  a  front  as  it  passes  a  layer.    We  make  two  simplifying  assumptions. 

1)  The  fractional  flow  function  is  the  same  in  both  layers  (i.e.  the  same  in  both  rock 
types). 

2)  The  flow  is  miscible. 
The  node  propagation  routine 


In  the  front  tracking  method  the  propagation  of  a  front  (in  the  hyperbolic  step  of  the 
sequential  scheme)  is  generally  accomplished  by  splitting  the  hyperbolic  operator  along  per- 
pendicular directions  in  a  local  coordinate  frame.  The  front  is  advanced  by  solving  a 
Riemann  problem  in  the  direction  normal  to  the  front  and  the  state  variable  (in  this  case  the 
saturation)  is  first  updated  according  to  the  normal  movement  of  the  front.  This  step  is  fol- 
lowed by  solving  the  hyperbolic  equation,  in  a  direction  locally  tangential  to  the  front.  This 
is  done  by  a  (one  dimensional)  finite  difference  scheme,  thus  updating  the  state  variable  at 
the  points  which  define  the  front,  and  accounting  for  flow  tangential  to  the  front.  Near  a 
point  of  interaction  of  fronts  (a  node)  these  steps  cannot  be  taken  in  the  usual  way  since  the 
problem  is  inherently  two  dimensional  and  splitting  into  local  normal  and  tangential  direc- 
tions is  no  longer  well  defined. 

An  algorithm  for  the  propagation  of  the  front  m  the  vicinity  of  a  layer  has  been 
developed  and  is  given  below.  Consider  a  point  on  the  front  that  is  to  be  propagated  past  the 
rock  discontinuity  (i.e.  from  the  upstream  side  (layer)  to  the  downstream  side  (layer)  of  the 
discontinuity).  The  algorithm  divides  the  propagation  into  two  parts.  First  the  point  is  pro- 
pagated using  the  component  of  the  velocity  normal  to  the  front  until  it  reaches  the  discon- 
tinuity. Then  an  angle  of  deflection,  velocity  and  the  new  normal  direction  for  the  propaga- 
tion into  the  downstream  layer  is  computed  using  only  the  information  in  the  upstream  layer. 
The  point  is  then  propagated  for  the  remaining  time  of  the  timestep  using  the  new  propaga- 
tion direction  and  velocity.  The  deflection  of  the  front  direction  at  the  discontinuity  is  analo- 
gous to  Fermat's  principle.  We  re-iterate  the  obvious,  namely  that  the  point  in  question  does 
not  represent  a  moving  fluid  particle  but  a  point  on  a  shock  surface  along  which  there  is 
tangential  slip. 

2.2.    Cross  flow 

As  an  application  of  the  algorithm  discussed  above,  we  study  a  case  of  interest  which 
occurs  when  the  flow  is  mainly  parallel  to  the  layer.  Consider  a  miscible  flow  initialized  as  in 
Fig.  2.1a  .  Assuming  that  the  left  layer  is  of  higher  permeability  than  the  right,  and  assum- 
ing no  cross  flow  between  the  layers,  at  a  later  time  the  solution  is  shown  in  Fig.  2.1b  . 
Cross  flow  between  the  layers  can  be  taken  into  account  qualitatively  by  noting  that  the  pres- 
sure distribution  in  the  no-cross  flow  solution  is  piecewise  linear  in  each  layer.  For  the 
unstable  flow  regime  (m  >  1)  these  pressure  distributions  are  shown  in  Fig.  2.1c  .  Thus  for 
y  <  c  the  pressure  distributions  would  favor  cross  flow  into  the  high  permeability  layer  and 
for  y  >  c  the  cross  flow  would  be  into  the  low  permeability  layer.  (See  [29].)  Under  the 
approximation  that  in  the  cross  flow  case  the  pressure  distributions  can  be  accurately 
represented  by  Fig.  2.1c,  it  can  be  shown  that 

c  —  a  1 


^-'         l  +  m(-^-   1) 
a 

where  L  is  the  length  of  the  computational  region  in  the  y  direction. 


(2.1) 


Fig.  2.2  shows  the  numerical  solution  of  this  problem  using  the  algorithm  developed 
above.  For  our  choice  of  parameters  the  ratio  (2.1)  is  about  0.1  .  The  direction  of  flow 
agrees  with  above  approximation;  over  the  major  portion  of  vertical  part  of  the  front  the 
cross  flow  is  toward  the  slow  region.  Therefore  one  expects  the  cross-over  point  C  to  be 
much  closer  to  the  slower  part  of  the  front  as  shown  in  Fig.  2.2  .  One  might  expect  that  the 
portion  of  the  front  between  A  and  C  should  move  into  the  fast  region,  as  shown  in  Fig. 
2. Id,  but  this  was  not  observed  numerically.  Even  when  the  front  was  initialized  as  in  Fig. 
2. Id,  the  front  did  not  persist  in  this  configuration. 

The  finger  formation  is  an  indication  of  a  singularity  at  the  node.  Initially  the  singular- 
ity becomes  stronger  as  the  finger  becomes  sharper  at  B,  thereby  accelerating  the  runaway 
behavior.  At  A,  the  reverse  happens;  as  the  angle  of  the  front  with  the  discontinuity  moves 
away  from  the  normal,  the  velocity  singularity  becomes  weaker. 


3.    Porosity  Variation 

3.1.    Introduction 

One  of  the  main  objectives  in  oil  recovery  is  to  suppress  the  fingering  and  channeling 
instabilities  which  are  initiated  by  small  and  large  scale  disturbances  through  the  nonunifor- 
mities  in  the  medium.  The  nonuniformities  that  we  have  in  mind  are  the  variations  in  the  per- 
meability and  the  porosity  of  the  medium.  In  [3]  the  fingering  problem  associated  with 
heterogeneity  in  the  permeability  field  was  studied.  See  also  [20]  for  a  discussion  including 
the  effects  of  capillary  pressure.  A  partial  remedy  to  this  effect  through  the  use  of  polymer 
flooding  was  analyzed  in  [4].  In  this  section  we  address  the  fingering  problem  associated 
with  variable  porosity. 

To  gain  an  insight  into  the  effect  of  porosity  variation,  consider  (1.1)  -  (1.3).  The 
effect  of  porosity  can  be  qualitatively  understood  as  follows.  If  the  porosity  were  constant,  it 
can  be  removed  from  (1.1)  (leaving  (1.2)  and  (1.3)  unchanged),  by  redefining  the  time  vari- 
able, t  -  t  =  t  /  <J).  As  a  consequence,  the  speed  of  any  wave  that  would  appear  in  the  solu- 
tion to  the  hyperbolic  equation  (1.1)  would  be  modified  by  a  constant  factor  inversely  pro- 
portional to  the  porosity.  In  the  case  of  variable  porosity  the  above  argument  can  be  applied 
locally,  to  leading  order  approximation,  thus  implying  a  spatial  variation  of  wave  speeds.  As 
in  the  case  of  nonuniform  permeability,  this  variation  can  act  to  produce  fingering. 

V 

If  we  introduce  a  new  velocity  field,   v=  — (x),    (1.1)  -  (1.3)  can  be  rewritten  as 

St  +  V  ■  (V  f(s))  =  f(s)  V  •  V  ,  (3.1) 

V  =  -  .^X  .  VP  ,  (3.2) 

<P 

V  •  KX  •  VP  =  0  .  (3.3) 
From  (3.2)  we  see  that  the  effective  rock  permeability  for  the  velocity  field  v  is  K  /  <|).    It  is 
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well  known  that  regions  of  local  maxima  in  the  permeability  field  serve  as  nuclei  for  finger 
growth  in  porous  media  flow  (see,  for  example,  [3|).  Thus  (3.2)  implies  that  regions  of  low 
porosity  will  have  an  effect  similar  to  that  of  high  permeability.  However,  since  the  porosity 
does  not  enter  into  the  elliptic  equation  (3.3),  its  effect  will  be  milder  than  variations  in  the 
permeability  field  itself.  In  addition  we  note  that  regions  of  high  reservoir  porosity  com- 
monly have  high  permeability.  Thus  porosity  variations  will  normally  provide  a  partial  offset 
to  the  fingering  tendency  of  permeability  variations. 

To  explore  further  the  effects  of  permeability,  we  consider  the  equations  (1.1)-(1.3)  in 
one  spatial  dimension.  This  is  particularly  appropriate  when  one  thinks  of  the  movement  of 
discontinuities  of  the  hyperbolic  solution  as  a  two  step  procedure  (as  in  the  algorithm  used  in 
the  front  tracking  method),  namely  the  propagation  of  the  discontinuity  in  a  locally  normal 
direction  followed  by  a  step  in  which  the  tangential  slip  of  the  fluid  along  the  interface  is 
accounted  for.  In  one  space  dimension,  the  seepage  velocity  v  is  constant  as  seen  from  ellip- 
tic equation  (1.3)  and  setting  it  to  unity,  without  any  loss  of  generality,  reduces  (1.1)  to 

(<i)s),  +  f(s),  =  0  .  (3.4) 

Expressing  (3.4)  in  terms  of  the  conserved  quantity  (t)S, 

(4)5):  +   -^(4)S)x=   --^^4),  .  (3.5) 


reveals  that,  along  the  characteristic  lines 

dx  fs 


dt  <i)(\) 


(3.6) 


(|)S  is  not  constant  but  changes  by 


d((l).s) 


f. 


Indeed,  the  effect  of  the  source  term  is  to  force  the  saturation,  s,  to  be  constant  along  charac- 
teristics, as  can  be  seen  by  expressing  (3.4)  in  the  non-conservative  form 

fs 
s,  +  -— -s,  =  0  .  (3.8) 

'        <|)(x) 

Clearly,  the  variable  porosity  affects  the  curvature  of  the  characteristics  in  space  time  (  equa- 
tion (3.6))  and  the  "time  to  shock  formation "  when  starting  with  smooth  data.  It  is  also  seen 
from  (3.6)  that  the  characteristics  will  not  be  smooth  at  a  point  of  discontinuity  in  (i>. 

In  [22],  source  terms  such  as  found  in  (3.5)  are  seen  to  lead  to  additional  standing 
waves  in  the  solution  of  the  one-dimensional  Riemann  problem  associated  with  a  hyperbolic 
equation.  Such  waves  do  not  arise  in  the  present  case  due  to  the  special  form  of  the  source 
terms  and,  in  particular,  the  constancy  of  s  along  characteristics. 
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3.2.  Code  modification  for  inclusion  of  permeability  effects 

As  the  variable  porosity  enters  into  the  system  only  through  the  hyperbolic  equation 
(3.1),  the  adaptation  of  the  front  tracking  method  to  the  case  of  variable  porosity  requires 
incorporating  its  effect  in  the  solution  of  the  hyperbolic  equation  only.  Before  we  discuss  the 
incorporation  of  the  porosity,  it  will  be  helpful  to  describe  briefly  the  basic  ideas  behind  solv- 
ing the  hyperbolic  equation  in  our  front  tracking  method.  At  any  fixed  time,  the  (bounded) 
spatial  domain  of  the  computation  consists  of  a  number  of  regions  in  which  the  solution  is 
smooth.  These  regions  are  separated  by  discontinuities  across  which  the  saturation  is  discon- 
tinuous. The  numerical  algorithm  to  advance  the  solution  of  the  hyperbolic  equation  (3.1) 
from  time  t  to  t  +  dt  is  done  by  a  spatial  splitting  of  the  hyperbolic  operator,  solving 
separately  for  the  propagation  of  the  discontinuities  (the  "front")  which  includes  the  solution 
for  the  saturation  immediately  on  each  side  of  the  discontinuities,  and  for  the  solution  in  the 
smooth  "interior"  regions. 

The  position  and  shape  of  each  discontinuity  is  resolved  by  a  finite  number  of  points. 
Each  point  of  a  discontinuity  is  advanced  in  a  direction  (locally)  normal  to  the  discontinuity 
by  solving  the  Riemann  problem  associated  with  the  normal  component  of  (3.1).  The 
Riemann  problem  solution  gives  both  the  saturation  information  immediately  ahead  and 
behind  the  propagated  point,  ar^d  its  speed  of  propagation.  The  solution  for  these  saturations 
is  dictated  solely  by  the  flux  function,  f(s).  Since  this  flux  function  does  not  depend  on  the 
porosity,  the  states  across  the  discontinuity  are  unaffected  by  the  variation  in  porosity.  How- 
ever, as  shown  above,  the  speed  of  propagation  is  proportional  to  the  inverse  of  the  local 
porosity.  Inclusion  of  porosity  effects  in  this  step  thus  consists  of  incorporating  the  porosity 
in  the  speed  of  each  front  point. 

The  tangential  flow  of  fluid  along  each  discontinuity  is  incorporated  through  the  solu- 
tion of  the  tangential  component  of  the  equation  (3.1).  This  is  accomplished  by  employing 
standard  one  dimensional  finite  difference  methods  (i.e.  upwind)  using  the  stencil  determined 
by  the  representative  points.  This  finite  difference  solution  is  obtained  separately  for  each 
side  of  the  discontinuity,  and  these  solutions  are  easily  modified  to  take  variable  porosity  into 
account. 

The  solution  of  the  two  dimensional  hyperbolic  equation  in  the  smooth  regions  uses  spa- 
tial \-y  operator  splitting  and  standard  one-dimensional  finite  difference  schemes  which,  as 
just  mentioned,  are  easily  modified  to  take  into  account  the  variation  of  porosity. 

3.3.  Test  problems 

We  compare  the  effects  of  a  variation  in  the  porosity  and  in  the  rock  permeability  fields 
on  fingering  in  an  areal,  quarter  five-spot  flood.  Let  x(\,  y)  be  a  gaussian  random  variable 
of  mean  one,  and  standard  variation  ct.  Let  (})  and  K  be  constant  values  for  rock  porosity 
and  permeability  respectively.  Then  4*  x(x,  y)  is  a  gaussian  random  field  for  the  porosity, 
with  mean  (j)  and  variation  ct  4).  A  similar  statement  holds  for  K  x(x,  y)-  Using  a  particu- 
lar choice  for  the  random  variable  that  allows  a  specification  of  a  given  length  scale  for  the 


-9- 

variation  (see  [20]),  we  consider  two  calculations  for  the  unstable  flow  regime  of  immiscible 
displacement;  the  first  with  a  variation  in  the  porosity  field  given  by  6  x^x.  y)  with  K  fixed, 
the  second  with  ^  constant  and  the  variation  in  the  permeability  field  given  by  K  x(\,  y). 
Fig.  3.1  shows  the  growth  of  fingers  due  to  the  variation  in  porosity.  The  local  maxima  and 
minima  in  the  porosity  field  are  shown  respectively  by  -i-  and  -  signs.  The  fingering  can  be 
seen  to  be  initiated  in  regions  of  local  minima  of  the  porosity.  Fig.  3.2  shows  the  effect  of 
variable  permeability  with  fixed  porosity.  The  permeability  variation  causes  much  stronger 
fingering  than  the  porosity  variation. 


4.    Shearer's  Theorem 

A  mathematical  analysis  [16, 17,26,27]  of  the  equations  for  three  phase  (oil,  gas,  water) 
flow  has  revealed  a  very  complex  pattern  of  nonlinear  waves  and  wave  interactions.  The 
no-go  theorem  of  M.  Shearer  states  that  under  very  general  hypotheses,  the  "worst  case" 
complications  occur  generically.  This  requires  a  careful  examination  of  both  the  hypotheses 
and  the  conclusions.    In  this  section  we  give  a  preliminary  analysis  of  the  conclusions. 

An  elliptic  region  must  occur  generically  in  the  three  phase  flow  equations,  according  to 
Shearer's  theorem.  Since  the  Cauchy  problem  is  ill-posed  for  elliptic  equations  and  since  the 
Cauchy  problem  must  be  solved,  an  elliptic  region  could  be  regarded  as  a  deficiency  in  the 
equations.  Here  we  adopt  an  opposite  point  of  view,  and  argue  that  one  can  learn  to  live 
with  the  elliptic  regime  [7]. 

A  mathematical  analysis  of  the  Riemann  problem  for  related  equations  reveals  no 
pathology  or  nonphysical  behavior  in  the  solution  [15].  A  numerical  solution  of  three  phase 
flow  equations  reaches  the  same  conclusion  [1].  The  elliptic  region  is  a  bounded,  interior  sub- 
set of  the  state  space.  The  elliptic  instability  appears  to  manifest  itself  by  causing  the  solution 
to  exit  from  this  region,  whereupon  it  enters  the  (stable)  hyperbolic  region.  Thus  the  equa- 
tions could  be  viewed  as  predominantly  hyperbolic  with  a  non-infinitesimal,  nonlinear  stabili- 
zation of  their  infinitesimally  unstable  elliptic  region. 

In  two  and  three  space  dimensions,  we  expect  this  hyperbolic  stabilization  to  behave  in 
the  manner  of  a  phase  transition  Ihe  fluid  will  prefer  to  flow  in  hyperbolic  portions  of  state 
space.  If  an  elliptic  concentration  of  phases  were  somehow  initialized,  we  expect  the  solution 
would  segregate  itself  into  spatially  coherent  blobs  of  mixtures,  with  each  blob  located  within 
the  hyperbolic  portion  of  phase  space.  For  a  first  order  phase  transition  described  by  a  single 
order  parameter,  the  coherent  blobs  (pure  phases)  in  state  space  at  the  edge  of  the  mixed 
phase  region  are  joined  pairwise  by  tie  lines.  The  tie  lines  then  uniquely  sweep  out  the 
mixed  phase  region,  and  any  point  in  the  mixed  phase  region  decomposes  into  the  pure 
phases  at  the  end  of  the  tie  line  it  lies  upon. 

In  the  present  case,  one  would  look  for  tie  lines  by  finding  a  unique  solution  of  the 
Riemann  problem.  However  the  reasoning  is  circular,  as  the  Riemann  problem  is  probably 
not  unique  in  exactly  this  region.    In  this  case  the  nonuniqueness  is  resolved  by  appeal  to 
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more  fundamental  equations,  including  capillary  pressure  diffusion  equations,  and  finally  tie 
lines  (or  some  more  complex  solution  behavior)  will  be  determined  as  a  consequence.  It 
would  be  most  useful  to  derive  an  associated  order  parameter. 

Since  the  elliptic  region  has  little  effect  on  the  solution,  it  is  tempting  to  "eliminate"  it 
by  deforming  the  equations  within  a  fixed  topological  equivalence  class.  What  are  the  result- 
ing hyperbolized  model  equations?  Evidently  if  there  are  tie  lines,  they  should  each  be 
shrunk  to  a  point,  leaving  a  line  of  umbilic  points  where  the  two  hyperbolic  wave  speeds 
coincide.  In  this  case  the  hyperbolic  behavior  in  a  neighborhood  of  the  elliptic  region  will 
coincide  with  that  studied  for  a  polymer  flood  oil  reservoir  [19,28,18].  This  latter  observa- 
tion has  been  made  previously  by  B.  Keyfitz.  However  for  parameters  typical  of  real  reser- 
voirs, the  elliptic  region  is  very  small.  Thus  the  above  line  of  degenerate  hyperbolic  points 
can,  in  a  further  approximation,  be  shrink  to  a  point.  Doing  this  leads  an  isolated  point  of 
degeneracy,  as  in  the  models  studied  by  Marchesin  and  coauthors. 


5.    Conclusions 

The  front  tracking  code  developed  by  the  authors  and  coworkers  has  been  subjected  to  a 
series  of  tests  in  the  petroleum  reservoir  application.  These  tests  are  continuing  and  are 
becoming  increasingly  representative  of  realistic  engineering  practice.  Fundamental  progress 
in  numerical  algorithms  (e.g.  the  front  untangling  and  bifurcation  algorithms  [8]  )  and  in 
mathematical  theory  (e.g.  the  solution  of  Riemann  problems  in  one  and  two  space  dimen- 
sions) was  necessary  for  the  success  of  these  tests. 
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Fig.  2.1  (a)  The  initial  setup  for  a  miscible  flow  run  through  two  rock  layers  separated 
by  a  sharp  boundary.  The  initial  oil-water  bank  is  horizontal,  the  higher  permeability 
layer  is  on  the  left,  (b)  The  (one  dimensional)  solution  assuming  no  flow  between  the 
layers,  (c)  The  pressure  solutions  for  (b).  The  solid  line  is  for  the  left  layer,  the  dashed 
line  for  the  right,  (d)  The  direction  of  flow  expected  from  (b)  and  (c)  if  cross  flow  is 
included. 
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Fig.  2.2      The  solution  for  the  problem  discussed  in  Fig.  2.1  computed  using  the  Front 
Tracking  Method. The  initialization  is  same  as  in  Fig  2.1(a). 
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Fig.  3.1     Growth  of  fingers  for  immiscible  displacement  in  a  heterogeneous  reservoir: 
the  temporal  evolution  of  the  oil-water  discontinuity.    The  +  (-)  signs  correspond  to 
local  maxima  (minima)  of  the  porosity  field.    The  viscosity  ratio  of  the  water  to  oil  is 
1  :  10. 
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Fig.  3.2    Growth  of  fingers  for  immiscible  displacement  in  a  heterogeneous  reservoir: 
the  temporal  evolution  of  the  oil-water  discontinuity.    The  +  (-)  signs  correspond  to 
local  maxima  (minima)  of  the  permeability  field.  The  viscosity  ratio  of  the  water  to  oil  is 
1  :10. 
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